Stability under radiation reaction of circular equatorial orbits around Kerr black holes 
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We examine the evolution, under gravitational radiation reaction, of slightly eccentric equatorial 
orbits of point particles around Kerr black holes. Our method involves numerical integration of 
the Sasaki-Nakamura equation. It is discovered that such orbits decrease in eccentricity throughout 
most of the inspiral, until shortly before the innermost stable circular orbit (ISCO), when a critical 
radius rcrit is reached beyond which the inspiralling orbits increase in eccentricity. It is shown that 
the number of orbits remaining in this last (eccentricity increasing) phase of the inspiral is an order 
of magnitude less for prograde orbits around rapidly spinning black holes than for retrograde orbits. 
In the extreme limit of a Kerr black hole with spin parameter a — 1, this critical radius descends 
into the "throat" of the black hole. 



I. INTRODUCTION 



. Gravitational waves emitted by solar-mass-size compact bodies orbiting massive (IO^Mq and greater) black holes 
0^ ' (and spiralling towards them as they lose energy and angular momentum to the emitted radiation) are a favoured 
source for gravitational wave detectors sensitive to low-frequency radiation, such as proposed space-based detectors 
like the Laser Interferometer Space Antenna (LISA) [Q . Systems of this type lend themselves to theoretical analysis 
" ' via perturbation theory, because of the extreme mass ratio between the two bodies. In recent years, the Teukolsky 
. perturbation formalism for black holes has been employed successfully to describe orbital decay of small bodies orbiting 
' a large Schwarzschild (i.e. non-rotating) black hole using a mixture of analytic and numerical results (for an 

T— I . informative overview of much of the purely analytic work to date, see Mino et al. ||]). One result of this work has 
' been to modify the long-standing result that, under radiation reaction, orbits tend to become more circular as they 
C3 , slowly decay. In fact, inside a critical radius, which is Tcrit = 6.6792M for nearly circular orbits in the Schwarzschild 
geometry, non-circular orbits tend to become more, rather than less, eccentric [Q. Although a precisely circular orbit 
would remain circular inside the critical radius, its circularity is no longer stable to small perturbations away from 
O precise circularity as the orbital decay continues. 

^ Despite their intrinsic interest, these results may prove of limited usefulness for any future low-frequency gravita- 

?H . tional wave detectors, since there is no reason to expect that large black holes should typically have no spin at all. 

Just the opposite (that they should exhibit strong rotation) is perhaps more to be expected 0. Therefore it is of 
J> , great interest to extend this type of analysis to the case of rotating, or Kerr black holes. This presents no difficulty 
for the Teukolsky perturbation formalism itself, which was developed for the Kerr metric, but a problem does arise 
in dealing with an additional constant of the motion which governs orbits around spinning black holes. Unlike the 
■ energy and angular momentum, whose flux can easily be determined from the waves far from the source, until very 
recently there was no clear understanding of how to calculate the amount of "Carter constant" carried away by the 
emitted radiation. In spite of this, it has been shown recently, for general orbits in Kerr, that circular orbits (defined 
as orbits of constant Boyer-Lindquist radius, and sometimes referred to as "quasi-circular") remain circular under 
radiation reaction |^-[ic|]. While progress continues in developing techniques for dealing with general orbits in Kerr 
pd]-|l3|, it now seems worthwhile to investigate the case of nearly-circular, equatorial orbits around rotating black 
holes |14|]. Equatorial orbits in the Kerr spacetime, like orbits in Schwarzschild, can be said to have zero "Carter con- 
stant" , which remains unchanged during orbital decay. Looking at these orbits can tell us if the behaviour previously 
observed for slightly-eccentric orbits in Schwarzschild is also seen in the Kerr metric for all values of the Kerr spin 
parameter < a < 1. 

It is shown in this paper that, for equatorial orbits, it is generally true that a critical radius, rcrit exists beyond 
which slightly eccentric orbits become less circular due to radiation reaction, and that this radius is encountered 
close to the radius of the innermost stable circular orbit (ISCO). This is best illustrated by examining the behavior 
of the parameter c = roe/er'o, where e is the orbital eccentricity, and ro the mean radius, and an overdot indicates 
differentiation by time. This parameter is negative for orbits evolving with increasing eccentricities, and positive for 
decreasing eccentricity. Near the ISCO one can show, as in Sec. 9 below, that c diverges towards negative infinity 
near the ISCO, for nearly all values of a. There is an apparent exception to this behaviour in the limiting case of 
a maximally rotating Kerr black hole with a = M. In that case, the horizon and the ISCO are both located at 
r = M in Boyer-Lindquist co-ordinates, although they are still separate in terms of proper radial distance. As one 
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approaches r = M, for the case of a prograde orbit around an extreme Kerr black hole, c is both postive and finite, 
approaching the limit of 3/2 at r = M . Not surprisingly therefore, for prograde orbits around black holes with very 
large a > .99M, the transition to eccentricity-increasing inspiral takes place only shortly before the onset of dynamical 
instability at the ISCO in terms of the Boyer-Lindquist radial co-ordinate. The number of orbits remaining at this 
point is an order of magnitude fewer for such cases than it is for retrograde orbits in the same geometry. 

Since the radius of the ISCO is much smaller for prograde than for retrograde orbits (with a = M, nsco — M for 
prograde orbits and nsco = 9M for retrograde orbits), the critical radius is also much smaller for prograde orbits. 
These results demonstrate that the onset of "back reaction instability" for circular orbits precedes, and is intimately 
connected with, the onset of dynamical instability signified by the ISCO. It seems reasonable to conjecture that the 
alteration in the shape of the radial protential as the ISCO approaches, at which point the minimum of the effective 
potential vanishes, is reponsible for the gain in eccentricity. 

The organisation of the paper is as follows. In section 2, the orbital equations for geodesic motion (i.e. not 
including radiation reaction) are solved analytically for slightly eccentric, equatorial orbits. In section 3 the Tuekolsky 
perturbation formalism is described, and section 4 shows how to calculate the fluxes of energy and angular momentum 
carried away from the system using this formalism. In section 5 the Sasaki-Nakamura equation, which is actually 
solved rather than the Tuekolsky radial equation for numerical reasons, is presented. In section 6 the Teukolsky 
source function is calculated for a perturbing particle following the orbits of section 2, and the results of both of these 
sections come together in section 7 to yield the rate of change of orbital eccentricity due to radiation damping. This 
orbital evolution is described under the assumption of adiabaticity (that the orbital evolution is much slower than the 
orbital period), which introduces constraints which are discussed in section 8. Finally, in section 9, the analytic and 
numerical results are presented, followed by a discussion of their significance in section 10. A guide to the essential 
points of the argument is given at the end of section 7. 



II. DESCRIPTION OF THE ORBIT 



Since the perturbation of the Kerr metric producing the gravitational waves is that of a small particle orbiting 
the black hole, it will be necessary to solve the orbital equations for a particle in orbit around a rotating black hole. 
We require expressions for r(i), (f)(t) and 6{t) to describe the orbit in Boyer-Lindquist co-ordinates. Since we restrict 
ourselves to equatorial orbits, the solution for the 6 motion is trivial, 9 = n/2 is a constant throughout. The equatorial 
orbital equations for a particle in the Kerr spacetime, in these co-ordinates (leaving aside the trivial dO/dr = 0), are 
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fiE'^dr/dr = [{E{r^ + a^) - aL,f - A{^?r^ + (L, - aEf)]^ = VR (1) 
fiY,'^d(j)/dT ^-{aE- L,/ sin^ 9) + {a/A){E{r^ + a^) - aL^) = $ (2) 

(2 I 2^ 

nT,'^dt/dT = -a{aE sin^ 9 - L^) + ^ — ^ — -{Eir"^ + a^) - aL^) = T, (3) 

where t is proper time, Yj = + cos^ 9, IS. = — 2Mr + a^, and the black hole's spin parameter a is defined for 
convenience as a = J - L/M , with J the spin angular momentum vector of the black hole, and L a unit vector pointing 
in the direction of the particle's orbital angular momentum vector. For prograde orbits (in which the particle orbits 
in the same sense as the black hole's spin) a is positive, and for retrograde orbits (in which the particle rotates in the 
opposite sense to the hole), a is negative. Recall that we restrict attention to equatorial orbits only, so that J and L 
are either parallel or anti-parallel. It is a condition of the perturbation scheme that fi/M ^ 1, where AI is the mass 
of the central black hole and /i the mass of the orbiting particle. Finally, E and are the particle's orbital energy 
and angular momentum, respectively. 

We now consider slightly eccentric orbits, and define a mean radius tq, so that d{R/r'^) / dr\r=ro = 0. The eccentricity 
e is defined so that R{r = ro(l -I- e)) = 0. These definitions are chosen so that as e ^ 0, rp reduces to the constant 
radius of a circular orbit, and so that e corresponds, when e <C 1 and in the appropriate limits, to definitions of the 
eccentricity of an orbit commonly used in the Schwarzschild geometry and in Newtonian mechanics . These defining 
equations for tq and e permit us to write the orbital energy and angular momentum in terms of these two quantities. 
Since we assume throughout that e is a small quantity, it is convenient to expand E and in terms of it, 

E{ro,e) = Eoiro) + eSi(ro) + e^E2{ro) + e^E^iro) + 0{e^) (4) 
L.{ro,e) = io(ro) + eii(ro) + r^L2{ro) + e^^L^iro) + 0{e^). (5) 

Using our two equations in vq and e, it is easy to show that 
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1 - 2t;^ + qv^ 

Ea = ^ (6) 

(1 - 3v^ + 2qv3)i 

Ei=0 (7) 

w^fl - + qv^ + g2x;4)(i _ gu^ + 8qv^ - 3q'^v*) 

E2 = u ^ (8) 

2(1 - 3i;2 + 2qi;3)i(l - 2i;2 + g2^4) 

w2(l-3w2+gu3_^52^4)(-]^_7^2_^;^Q^^3_4^2^4) 

-C/3 = — u o (9) 

(1 - 3t;2 + 2(7w3)i(l - 2i;2 + g2„4) 

Lo ^ /o-(l-2g-^+-^^f) (10) 
(l-3t;2 + 2qw3)5 

Li=0 (11) 

^ _ qrov^{q-3v + qv^ +q^v'^){l-&v'^ + Sqv'^ -'iq^v'^) 

2(1 - 3i;2 + 2gw3)i(l - 2i;2 + q^v^) 

qrcv^iq-iv + qv"^ +q^v-^){l-7v'^ + lQqv''^ -Aq^v^) 
L3 — — /i 5 . (13) 

(1 - 3i;2 + 2qv3)i(l _ 2w2 + g2„4) 



Here v = ^JM/tq and q — a/M. These results, up to order are given in Ref. |l4|] . 

We wish to write the change in the eccentricity brought about by the loss of orbital angular momentum and energy, 
in terms of the rates of loss of those two quantities. Since we have E and Lz as functions of and e, we use the 
chain rule for differentiation to write 

dE BE 

E = -dEcw/dt = — e + —fo (14) 
oe oro 

dL dL 

Lz = -dLcw/dt = -^e + -^fo, (15) 

where dEaw / dt and dLaw / dt are the total energy and angular momentum carried towards infinity and the black hole 
horizon per unit time by the gravitational waves, averaged over several wavelengths. We will write these quantities 
also in terms of e and tq, 

dEaw ^ , , , ^,3 



dt 



^ Eo + eEi+e'E2 + 0{e-') (16) 



dt 

As we shall find later, Ei — Li — 0. Eliminating ro from Eq. (ffsh, we derive 



^^=Lo + eL,+e^L2 + Oie'). (17) 



dEgw , dLgw -p,,.i,dE ^L^ , 



where / = d/drQ. 

Substiting Eqs. ( pT| ) and (^ into Eq. ([l^), we find, keeping terms up to order e2, 

-L'.iEo - #io) - e^L',iE2 - ^U) - e^L'^iEo - ^U) 
■ _ ^ ^ ^2 

2e{E2L'^ ~ L2E'^) 

Now, from Eqs. (|l3|), we see that 

E'r, VM 



(19) 



7T = 3 '"^ ^ - f^, (20) 



^0 + aVM 

where f2 is the angular frequency of a circular orbit of radius rg. It follows from an interesting (and quite general 
1^ ) characteristic of circular orbits, and will be shown later in this case that, the circular (i.e. zeroth order in the 
eccentricity) rates of loss of energy and angular momentum are related by 

^0 = nU. (21) 
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Therefore 



where 



and 



where 



^ie^-eJ{v)[g{v)Eo + E2-nL2], (22) 



_ /i _ (1 + qv^)il - 2v^ + q^v^){l - + 2qv^)2 



niv)^^-^^ ^-^ (24) 

^ L'q E'q 2(l + gv3)(l-6w2 + 8qw3_3g2^4)(i_2i;2 + g2„4)2' V ) 



g{v) = 2 - 27v^ + 12v^ - + 38qv^ - Uq^v"^ ~ lUqv^ + 86q^v^ 
+ iq^v'^ + 72qv'^ ~ 12q^v^ - 36q^v^ - 23q^v^" + SOg^i;" 

- 9g%i2 (25) 

Since e is proportional to e in this equation, it is plain that a precisely circular orbit (one with e = 0), will remain 
circular under radiation reaction, provided that we can indeed show that Eq = ^ILq and Ei = Li — 0. It is also plain 
that the question of the stability of an orbit's circularity will be determined by the sign of Eq. (^^, which requires 
us to calculate the loss of orbital energy and angular momentum up to second order in e. 

Similarly we can solve for fg, the rate of change of the orbital radius, which tells us that to leading order tq = 
—Eq/Eq, which implies that 



2(1 - 3t;2 + 29t;3)3/2 
(1 - 6v^+8qv^ - 3(j2 ^4 



f^ro/ro = -— — TTir-^ri:^ — (^6) 



With this in hand it is possible to proceed to the solution of the geodesic equations [Eqs. (^]. We expand r{t) 
about the mean radius tq in terms of the small eccentricity e, so that 

r{t) = ro[l + eri{t) + e^r2{t) + 0(6^)]. (27) 

Making use of the expansions of i?, and r{t) in terms of e, we expand out the equation (dr/dt)'^ = R/T^, and 
collect terms of order and (note that the term in r{t) does not contribute until 0{e'^) in R/T'^), giving us two 
differential equations, 

{dn/dtf = nl{l~rl), (28) 

where we define a radial frequency, 

nr^fl{l~6v^+8qv^~3q^v'^)i (29) 



and 



where 



ilr^ + ^^1^2 = fi{v) + hivVi + h{v)rl (30) 



1 - 7z;2 + IQqv^ ~ Aq^v^ 
^'^""'^ l-6v^ + 8qv^-3q^v^ 
, . ^ 2t;^(l-2gt;3 + g2^4) 

il + qv^){l-2v^ + qH^) ^'^ > 

■'^^ ^ {\ + qv^){\-2v'^ + q^v^){l-Qv-^ + 8qv^ -3q^v^y ^ ' 



and 
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- IQqv'^ + Iq^v^ + 2Aq^v^ - Aq'^v^ - 27q^v^ + 16qV° 
-4gV\ 



Integrating these equations in order, we find, 



r2{t) = -h{v){l - cosiSlrt)) + -h{v){l - cos(2f],t)). 



(34) 

(35) 
(36) 



It remains to solve for the (/)-motion. Again we expand out the geodesic equation d(/)/dt = ^/T, integration of which 
yields 



where 



and 



p{v) 



(l){t) = n^t - ep{v) sin(f^^t) + 0(6^), 
2(1 - 3v^ + 2qv^) 

[(1 + qv^){l - 2v^ + q^V^){\ - 6^2 + SgV^ _ 3g2^4)l/2] 



(37) 



(38) 



3(1 - 111)2 + 24u'* + liqv^ - 4g2v4 _ 46qi,5 + 2hq^v^ + q^v^ - iq^v^) ^ 



2(1 + gw3)(l - 2z;2 + g2„4)(]^ - 6i;2 + ; 
0[1 - A17e2 + 0(e3)] 



3g2i;4) 



(39) 



is the azimuthal angular frequency. The 0(e2) part of (/)(t) which is proportional to sin(Ori) is not given, as neither 
it nor the 0(e2) part of r{€) contribute to the final result for e, for reasons which will become clear later. Only the 
0(e2) part of fi^ (i.e. Afi) is required, although it is necessary to know r2(t) to derive ASl. 



III. THE TEUKOLSKY FORMALISM 



We employ a scheme previously used in the Schwarzschild case to investigate the evolution of slightly eccentric 
orbits under radiation reaction m. This scheme is based on the Teukolsky formalism for perturbations of the Kerr 
metric. In this formalism one can decompose the Weyl scalar ^/)4 (which describes gravitational wave fluxes near 
infinity for such a system) as follows, 

= (.-.acos.)4 y_ E ^->^(^)-^r™(^)e'"^e--d., (40) 

where -28^^ is the spheroidal harmonic function of spin weight s ~ —2. The normalization used here for these 
fmictions is \-2S^j'^^{9)\^ sinOdO ~ 1/2tt. The radial function Rimuj{r) obeys the Teukolsky equation, 

A^^(^^) -y(r)i?,™.(r) =T,_(r), (41) 
where Timui is the Teukolsky source function, to be evaluated below. The Teukolsky potential is defined by 

V{r) = — ^- — + d,iLur + X, (42) 

where K = {r"^ + 0^)0) — ma and A is the eigenvalue associated with the appropriate spheroidal harmonic _2'S'°m- 

We can define two solutions to the homogeneous Teukolsky equation, Rf^^ir) and -R;^„(f), with the following 
boundary conditions, 
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and 



i?;^^^AV'='-*,as r^r+ (43) 
RLu.-r'BZle'-^' + -BZ^e-^-'',^s r ->(X) (44) 



RZu.-r'(^-'^'^\^s r^oo, (46) 



where A; = a; — ma/(2Mr+), r+ = M + V Af ^ — is the radius of the black hole horizon, and r* , the tortoise 
co-ordinate, is defined as 

r*=r+ ^In In , (47) 

r I - r_ 2AI r. - r_ 2M ' ^ ' 



where r_ = Af — v Af ^ — . 

From Ref. [0, the solution of the Teukolsky equation (solved via a retarded Green's function) is 

Rim^ir) = RZMZ"{r) + Rg,Jr)Z°-{r), (48) 

where 

r+ ^ 



Z"{r) = ^ / 



Z~(r) = \ / (50) 



and 

For convenience, we will write Zf^^^ ~ Z^{r oo) and Zj^^^ — Z°°{r ^ r+), and therefore our two solutions can 
be written as 

Rin^Ar-^^)- Zl'^.ye'^^' (51) 

and 

RirnUr ^ r+) ^ ZZ^A\-^'^^' . (52) 



IV. ENERGY AND ANGULAR MOMENTUM FLUXES 



Towards infinity, the Weyl scalar can be related to the two fundamental polarizations of gravitational waves by 

ip4 ^ ^(h+ - i'hx)- (53) 

From this and Eq. (|40| ) above, we can determine the averaged energy and angular momentum fluxes at infinity, 
employing the Isaacson stress-energy tensor to define the energy flux in the wave , as 

dEg-w _ poo _ l^imfcP 



and 



lmk 



lmk 

where the amplitude coefficient is decomposed into a discrete set of frequencies based on the particle's orbital motion. 
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ZlL^^Y.^"^k5{i^-^'^)- (56) 
fe 

Energy and angular momentum are also lost by radiation through the horizon of the central black hole. Again, ip4^ 
completely describes the waves as r* —oo and, with Teukolsky and Press we find 

E" = Y^af\^ (57) 



and 



with an identical decompostion of as with Zf^^, and where 

_ 2'^wlk{P + 4£')(fc' + 16e')(2Mr+)5 



and e = - a?/AMr+ and 



(59) 



|Cp = [(A + 2f + Aaum - Aa^uj\\^ + iQaum - SGa^w^) 

+ (2A + 3)(96a2tj2 „ A^aum) + 144^2(^^2 _ ^2^^ ^-gQ^ 

The total rates of loss of energy and angular momentum by the system are + and + L^. 



V. THE SASAKI-NAKAMURA EQUATION 



The preceding section makes it clear that our chief task is to calculate the amplitudes Zf^"^ , and it is apparent 
from Eqs. (^) and (^) that this will entail solving the Teukolsky equation to find the amplitude of the in-going 
waves at infinity, B^'^^^ from Eq. (^^. Numerically this presents a problem, however, since the ingoing waves for this 



solution are completely swamped by the outoing waves at large radii [compare amplitudes of 



and Bll^^/r as 



r oo in Eq. (^)]. In the Schwarzschild case this problem is typically avoided by solving instead the Regge- Wheeler 
equation, and transforming its solution to that of the Teukolsky equation via the Chandrasekhar transformation [po[ . 
The virtue of this is that in the Regge- Wheeler formalism, with its short-range potential, the ingoing and outgoing 
waves near infinity have the same order of magnitude. 

In the Kerr case Sasaki and Nakamura have found an equation with the same useful properties as the Regge- Wheeler 
equation in Schwarzschild which, moreover, reduces to the latter equation when a —> |21|. The Sasaki-Nakamura 
equation is written as follows 



Fir] 



dr 



(61) 



The functions F{r) and U{r) are given in the appendix. The equivalents to our two solutions to the Teukolsky 
equation are 



Imuj 



X, 



H 

Imuj 



/tout iojr* 

e^'^''' , as 



A 



linuj 



'r+ 



oo 



(62) 
(63) 



and 



A;. 



IrriLU 



e , as 

jjQut ^ikr* 



(64) 
(65) 



The transformations between the quantities we require are 



7 



and 



oin _ ^ A in 

where Xhrmj = -^im^^/"^^^ + ^"^"^ "^Oi cij |5 and 77 are given in the appendix. 

VI. THE SOURCE TERM 

The Teukolsky source term is given by 



where 



Ti^u. = 4 J p-^p-\B'2+B'2*)e-'"''^+''^'-2SZdndt, 



B'2 - -ip«pL_i[p-4io(p-'p-^r„„)] 

^ /pA2i_i[p-4p2j+(p-2p-2A-i7;,,„)], 



and p = {r — ia cos 9)^^ , with p its complex conjugate. The operators Lg and J+ are defined 

TTi 

Ls = H — ^ — - — aui sin 9 + s cot 
sin W 

and 

J+ ^ dr + l-^. 

The tetrad components of the particle's energy momentum tensor can be written 



T — 


Clin c { 

smt' 


r{t))6{9 ~ 


7r/2)<5(0- 




T- — 




r{tm9- 


7r/2)5(0- 




T- - — 


sine ^ 


~r{t))5{9^ 


-7r/2)(5(0- 





where 



'-'nn — '-'nn ^ ^nn ^ '-'nn ^ ^-j- ' 



pi fdr ^2 



'-'mn — '-^mn ~r ^mn 



2v2S'^i sm i 



_ _^[jsin0(ai?--^]^ (78) 

Cmr-a^^[^sin6{aE^J^)f (79) 
zLr sm 

and t = dt/dr. 

Integrating by parts, and making use of the adjoint operator Ll = de — m/ sin9 + auj sin 9 + scot 6 — dg + f{9), 
which bears the foUowing useful relation to the operator Ls defined above: 

h{9)Ls [g{9)] sin 9d9 = - f g{9)L\_^ [h{e)] sin 9d9, (80) 
Jo 

with g{9) and h{9) arbitrary functions we find that 

)5(r-r(i)} 

-oo Jo 

+ {(^mril + A^iml)5{r - r{t))) 

+ {A,,,rn25{r - r{t))}^rrW " Ht)y^'-'""^d(t>dt. (81) 

The ^'s are all functions of r only, and in each case A — A^"^ + A^^^dr/dt) + A^^^dr/dt)^, where (writing -2<S'/'^ 
simply as 5" hereafter for simplicity) 

41 = -^Ci'y{rS,0e - 2iaSM + 2rf{TT/2)S.e 

-2iaf{Ti/2)S + rS{f{n/2f -2), (82) 

A^na = '^-^ct\r\S,e + /(7r/2)5)(z| + -J, (83) 



r 



A^^ - -r'C^r.S{-^{^),r - {^ f + f f ), (84) 
A^ni = ^r'C^iS,9 + f{n/2)S), (85) 

A^il, = -2r^Crn^Sii^ + -), (86) 
A r 

AfJm2 = ^T^'^CfnmS, (87) 

4(2) _ .(2) _ .(1) _ .(2) _ .(1) _ .(2) _ ^(1) _ J2) , . 

In every case the spheroidal harmonic function (S) and its derivatives are evaluated at 9 = -k 12. 

It is now easy to show, from Eqs. (^9|) and (^0|) and using integration by parts (keeping in mind that we are 
interested only in closed orbits, for which r+ < r < oo always holds strictly), that 

-1 /'OC /'OO /'27T 

yH.OO ^ ill tH,00, 



Zf:^ = / / / lL^(.r)5{r - r(tm^ - 0(t))d0dMr, (89) 

2^^^lmuj Jr+ J-oo Jo 

for which /,^~(r) = Cjr) + I^±{r)idr/dt) + ll^Jr){dr/dtf, where 

jT)H,oo 72 n^:C<3 

It is necessary to expand Z^^'^ in terms of the eccentricity e, keeping in mind that we wish, as shown in section 2 
above, to find E^ °° and L^'°^ to second order in e, and that each of these is proportional to However, it 

transpires that only terms up to order e in the integrand of Eq. ( |89| ) contribute to order in e, the rate of change of 
eccentricity derived from E^^°° and L^'°° . The reasons for this emerge as we proceed to expand Ij^^ (j), Sir — rit)) 
and 5(0 — 0(i)) in powers of e. 
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Employing the expansions of r(t) and (j)(t) derived above [Eqs. ( p7| ) and (|37|)], we can write the product of delta 
functions in Eq. ( |89| ) as a product of two Taylor expansions in the small parameter e, about the points r — tq and 

(j) - n^t. 

5{r - r{t))5{<j) - (t){t)) = 5{r - ra)5{(l) ~ fl^t) - ero cos n^ir - ro)5{(t) - fl^t) 

- ep{ro) sinnrtd' {(t) - n^t)S{r ~ ro) + 0{e^), (91) 

where the prime denotes differentiation with respect to the function's argument. We can integrate by parts in Eq. 
(^) to integrate terms containing derivatives of delta functions, and this will simply mean that 6'{(f> — fl,f,t) will be 
replaced by im5{4) — ^,(,1), since e"*""^ is the only other part of the integrand which depends on (j). Completing the 
(f> integration thus leaves us with the overall factor e'('^~™^'*)*, and some terms depending on cosQrt, sinQrt and, in 
the O(e^) part, on cos 2nrt and sin 2nrt. Following the time integration, then, we will have a series of delta functions 
of the type d{uj — mil^) [at all orders, except 0(e)], ^(ci; — mil^ ± fir) (at all orders from 0{e) up) and, in general, 
S{uj — mVlff, ± fcl^r) at 0{e^) and above. 

These delta functions, after integration over uj to derive 7/14 [Eq. (^)], produce terms representing energy and 
angular momentum emitted at the fundamental (circular motion) frequency Wm = mQ^, and at a series of discrete 
sidebands, w± — mfl^±flr and w±k = mfl^±kflr. The occurrence of these delta functions justifies the decomposition 
of Zf^'^ referred to earlier [Eq. ( p6| ) above]. 

It is, of course, \Zf^'^\'^ which is integrated in Eq. (|o|). Therefore, up to order e^, only those O(e^) terms in Zf^"^ 
which cross multiply with 0(6*^) terms will contribute. Since the frequency must be single valued for any given term, 
only the circular harmonic {w,n) term in O(e^) survives the Fourier transform which produces the Weyl scalar, all 
other terms being annihilated. The 0{e) terms in Z have no circular harmonic term, as mentioned before, so these 
terms only contribute to loss of energy and engular momentum at O(e^). 

As seen from E g. (2^ ) above, it is the difference E2 — ^1-^2 on which e actually depends at leading order. Eqs. 
(||),(|5|),(|^) and (Hfylhow that 

En - nin cx 1 , at order e" (92) 

which is zero to leading order if cj^ ~ utni ~ mfl^. This means not only that the O(e^) terms in Z^^ do not 
contribute at all to e below 0{e^), but also that E^ — flLo is also zero to leading order, as noted above [Eq. (pl|)]. In 
fact, since the eccentric correction to the azimuthal frequency 51^ is itself of O(e^), the circular losses of energy and 
angular momentum contribute to e at O(e^) to leading order, like the first order terms in Z. Therefore there is no 
loss of E and at 0(e), and so i^i = and Li = as claimed in section 2. 

This proves that a precisely circular equatorial orbit in Kerr will always remain circular under radiation reaction 
(as long as the adiabatic approximation still holds). Furthermore it means that to find the leading order correction 
to this condition for slightly eccentric orbits, and thus establish the stability of circularity, we need only examine the 
0(e) terms in Eq. (^0|), and can drop all O(e^) corrections to the motion, except for the Ail part of fl^. This also 
means, of course, that only contributions to the loss of energy and angular momentum from the first pair of sidebands 
{lu = Lu±) need be included with the circular harmonic (oj,,,.) in calculating e to leading order. 



VII. CALCULATION OF RATE OF CHANGE OF ECCENTRICITY 

As a final step before integration of Eq. ( ^9| ) , the function I^^]^ (r) must be expanded up to first order in e. It contains 
terms which depend on dr/dt which, by Eq. (^7|) above, is 0(e) at leading order, dr/dt = — erpfir sinfi^i + O(e^). 
Therefore we will write 

M = Cir) - el^l{r)r,nrsinnrt + 0{e% (93) 
Thus, doing a final integration by parts in the integral over r in Eq. (^9|), we find 

= - run,) - eB+J{u: - ^+) - eB-J{^ c._) + 0{e')], (94) 

Imuj 



where 

,t(0) 

± mp{ro)Cjro) T /;,l(ro)rof^.). (95) 



1 dT^"'> 
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The argument of the preceding section shows that, in order to calculate the quantity E2 — ^11^2 + A57ijo, we need 
only evaluate the co-efficients in Zim- Therefore, returning to Eq. (|2^), we have 

^ie/e^~Jiv)[T~h{v)EQ] (96) 

where 



r = £^2 - + AnEo (97) 



,3 



«r-^^«r (98) 



and 



with 



h{v) ^An- g{v) = ^ ^^3^^^ _ 2^;2 + g2„4)2(i _ 6^2 + 8^„3 _ 3^2„4) ' (^9) 



7^(z;) = 1 - 12^2 + 66w* - 108^^ + qv^ + 8q^v^ - 72qv^ - 20q^v^ 

-27g^t;". (100) 

As an aside, we take the opportunity to write the eccentricity in terms of quantities which can be deduced from the 
signal observed in a detector such as LISA. The complex wave amplitude at earth h{t) ~ — ih^ can be written as 

m = - ^ ^^f™.-25.r(^)e™^e--(--)d., (101) 

where r is the distance from the source to Earth and t ~ r* is retarded time. A glance at Eq. ( |94| ) tells us we can 
define, based on this equation, amplitudes for the main sideband with frequency ujm (call this amplitude hm) and 
for the various sidebands (let hi be the amplitude for the sideband of frequency (j-'+). To leading order hm will not 
depend on the eccentricity, whereas hi, the amplitude of the first sideband, will be linear in e. It is therefore easy to 
show that the eccentricity will be proportional to the ratio of the amplitudes of the first and the main sidebands (i.e. 
hi/hm)- In fact. 



hi 



(102) 



E ir^i?,+_,(ro)-25r:M^)e^™^e— ' 

In order to measure e as it evolves with the signal, the signal will have to be strong enough to permit not only 
measuring the size of the first sideband, but also some parameter extraction, so that a and M can be estimated. 



In summary, Eq. (96) is the equation which allows us to compute the change in eccentricity for an inspiralling orbit, 
and Eq. (|6|) defines the rate of inspiral. Eq. dsj), Eq. (|9|) and Eq. for T, and Eqs. (||) and (|^) with the 

0(e°) part of Eq. ( p4|) for Eq, allow us to express e in terms of the solution of the radial Tuekolsky equation 
and its derivatives, as well as the incoming wave amplitude Bl'^^. These quantities are in turn derived numerically 
by solving the Sasaki-Nakamura equation as described below in section 9, and employing the transformations given 



in Eqs. (pq), (67) and (pq). The important functions j{v), h{v) and Ail in Eq. (M) are all derived in solving the 



equations of geodesic motion for the orbiting body in section 2. 



VIII. ADIABATIC CONDITION 



The whole preceding argument depends on an adiabatic condition on the motion which says that the inspiral 
timescale ro/|ro| is much greater than the orbital period of the motion 27r/rir- The necessity for this condition is most 
noticeable in the approximation which describes the evolution of the particle's motion under back reaction as passing 
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through a series of geodesic orbits, each defined as if no back reaction were taking place during that orbit. Once the 
inspiral proceeds on a timescale which is about as short as the time to complete an orbit, this approximation loses all 
validity. Using Eq. (|6|), we find that the adiabatic condition can be written, 

/i (1 - + 8gu3 _ 3g2^4)3/2 ^ 

Just as the inspiral timescale must be greater than an orbital period, so too must the circularization timescale e/e. 



However, this quantity is almost invariably less than the inspiral timescale, so Eq. (103) is the key condition. For 
very large radii, in the Newtonian limit, {M/ ^)'^Eq ~ 32u^"/5 (for a discussion of this limit see Ref. (1) and the 
condition is simply ii/M <C (5/1287r)w~^, which is very much less restrictive than the linear perturbation condition 
/i/M ^ 1, upon which the Teukolsky formalism rests. Approaching the ISCO however, where the numerical results 
tell us that Eq remains finite and of the same order as its Newtonian value, we see that the adiabatic limit on ^./M 
is proportional to (1 — + Sqv^ — 3q^w^)^/^, which becomes vanishingly small as the ISCO nears. Therefore, near 
the ISCO the adiabatic condition supercedes the linear perturbation condition, as the leading constraint on /i/M. 
Only by imagining a test particle which has vanishingly small mass can we apply the results of our calculation all the 
way to the ISCO, but no doubt there exist real physical systems, with ii/M < 10~^ for instance, which are correctly 
described for almost all of the inspiral by this approximation (recalling that our calculations presume that the particle 
is a point mass as a further simplification). This issue will be discussed more quantitatively in Ref. [p3[. 



IX. RESULTS 

With the results of section 7, it only remains to calculate i?,^^, B}"^^^ [Eqs. (B|) and (|4|)] and _2S'J'^(7r/2) [Eq. 
(^)] numerically to find e/e. To find the solutions to the radial equation [Eq. (41)] one actually solves the Sasaki- 
Nakamura equation [Eq. (|6l|)] for X^[^ and A™^^ [Eqs. (|6^ ) and (|6^)]. These solutions are very smooth, apart from 
a singularity at the horizon r+ , and so Bulirsch-Stocr integration works very well in integrating them. The singularity 
is avoided by starting the integration from a point just outside the horizon (typically at + 10^^). The solutions are 
insensitive to variations by several orders of magnitude of this small increment. Richardson polynomial extrapolation 
is used to evaluate A^^^^ as r oo, since it can be expressed as the first term in a polynomial in 1/ujr defining 
the amplitude of the ingoing wave at large r in Eq. (^ [|j . This amplitude is evaluated for several endpoints of 
integration, doubling the endpoint radius at each trial, allowing the extrapolator to evaluate the limit of the amplitude 
as r ^ oo, which is 

The Spheroidal harmonic functions are calculated by expressing them as a linear combination of spherical harmonics 
of equal m, summed over all available values in I (truncating the series after 30 terms in practice). Substituting this 
series into the second-order ODE defining the spheroidal harmonics gives us a 5-term recurrence relation for the 
co-efficients of the expansion. This procedure, for the scalar case only, is found in |^6| . The recurrence relation for the 
expansion co-efficients can be solved using matrix eigenvalue routines which, like the Bulirsch-Stoer integrator and 
the polynomial extrapolator, are found in Ref. |^^ . The derivative of each spheroidal harmonic is also expressible as 
a combination of spherical harmonics of different spin- weight values by use of the edth operato r |B5| . Useful checks 
for the numerical results are found in the Schwarzschild limit, in j^] and in the circular limit, in [^7|. Analytically the 
results of sections 2 and 7 reduce to those of [^ in the Schwarzschild limit and those of section 2 to the results of p^ ] 
in the post-Newtonian limit. 

The accuracy of the numerical results is limited by several factors. The relative accuracies of the Bulirsch-Stoer 
integrator and the Richardson extrapolator can be increased easily, at some loss in computing speed. For these 
calculations they were set to 10~^ and 10~^ respectively. The solution of the eigenvalue problem has very good 
accuracy, but the approximation of the spheroidal harmonics as a combination of spherical harmonics begins to lose 
accuracy seriously when auj becomes much larger than order unity. However, this only occurs for very high (m > 20) 
harmonics of the motion for small radii, and these contributions are not required at the accuracy used here. The 
chief limit on accuracy is, in fact, the number of harmonics in I and m which are calculated. Invariably, for small 
eccentricity orbits, the leading order contribution is for / = 2, to = 2, and the significance of the contribution decreases 
sharply (but less so for small radii) with increasing I and m. A simple estimate, used in Ref. enables one to reliably 
estimate the inaccuracy involved in truncating the calculation at I = ^max- It tells us that, for a relative error (in 
estimates of the loss of energy and angular momentum) no greater than 77, with a mean orbital radius tq, then 
^max > logry/ log(M/ro) -I- 3. Taking all of these factors into account, we can generally estimate the accuracy of the 
numerical results at 10~^, and certainly the relative errors should be no greater than 10"'^ in the worst case. 

A useful parameter with which to investigate the orbital evolution is c, which represents a ratio of the inspiral 
timescale to the circularization timescale, or 
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ro de/dt 
e dro/dt 



(104) 



Again, c is positive when radiation reaction circularizes the orbit, and negative when it drives the orbit more eccentric. 
In order to see analytically the behaviour of c as the ISCO approaches, recall Eq. (pq) and write 



j{v)[T ^ h{v)Eo] 



(105) 



As ro 



I)] diverges, since rj^g^o - 

0. Since the numerical results show that F and Eq remain finite in all cases, it 



nscOi the radius of the innermost stable circular orbit, the function h{v) [Eq, 
3a2 



6Mrisco + Say/Mnsco 

is apparent that F (which is otherwise dominant), contributes negligibly near rigco- Therefore, making use of the 
expression for ro from Eq. (Ea), we find for ro near nscO: 



n 



4(1 - 2v'^ + qv^){l - + 2qv^){l - 6^2 + 8qv^ - 3q^v'^) ' 



(106) 



Again, 1 — 6v^ + 8qv^ — Bg^ti^ ^ Q as r ngco, so c diverges at the ISCO. However, its sign as this point approaches 



depends on the function Ti. [Eq. (100)], since the expressions in the denominator are all positive for r > rjsco. It is 
obvious that for large r, Ti. is always positive, but for small values of r, which can be achieved by prograde orbits 
around rapidly spinning black holes (a > .95M), Ti can become negative. However, it always becomes positive again 
before the ISCO, so that c ^ — oo at the ISCO, in all cases except one. 

The e xcep tional case is the extreme one of a — > M. A t this unique point, Ti., and all expressions in the denominator 
of Eq. ( |l06|) go to zero. Setting q = 1 in Eq. (106), and canceling factors of (u — 1) from both numerator and 
denominator, one finds that 



lim 



3/2, 



(107) 



which is both positive and finite, in contrast to the usual behaviour as the ISCO approaches. 

As Fig. 1 shows, the curves describing the critical radius and the ISCO do approach each other in terms of the 
Boyer-Lindquist radial coordinate as a ^ M, as our analysis of c might suggest. Therefore it is interesting to 
investigate the consequences of this for massive particles inspiraling around near extreme Kerr black holes. A useful 
measure here is the number of orbits left in the inspiral once the particle reaches the critical radius, that is, the 
number of orbits it will take the particle to reach the ISCO. Defining tc as the inspiral time between rcrit and nscOj 
and referring to Eq. (Eq) for the rate of inspiral, we have 



tr 



V^{1 - + 8qv^ ~ 3q2u4) 



2Eo{l- 3w2 + 29^3) 



(108) 



To a rough approximation, we can take Eq as constant in this region, and therefore 



JJ_ 

Eq 



1 



1 



v/l-3.2 + 2,.3(__ + __^^^^^^^^^ 



Approximately, the number of orbits left in this time will be 



tc 



1 



T 2tt M 2ttEq 1 + qvin 



^-^^^ + ^^^^'-2 ■ 2(l-3.2 + 2,.3) 



(109) 



(110) 



Note that Eq oc (/i/M)^, so that Nc is inversely proportional to fj,/M. In the test particle limit fJ,/M 0, Nc oo. 

For a = -.9M, we find that « .035M/yLt, while for a = .99M, Nc w .0025M//^. Note that the rate of energy 
loss is similar in these two cases (retrograde orbits radiate more energy for an orbit of given radius than do p rogr ade 
orbits), but the distance between rcrit and ngco is much smaller in the latter case. The condition of Eq. (103) at 
the critical radius for a — .99M is /i/M ^ .01, so these estimates are still applicable to systems with extreme mass 
ratios, such as compact solar-mass-size objects spiralling into rapidly rotating supermassive black holes. For such a 
system, a prograde orbit spends an order of magnitude or more fewer orbits in the eccentricity increasing phase than 
does a retrograde orbit. Furthermore, the orbital periods for these two cases (a prograde orbit with ro ~ 1.5Af, and a 
retrograde orbit with ro 9.5Af) are also very different, with the period of the retrograde orbit an order of magnitude 
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longer. The retrograde orbit therefore spends a factor of hundreds more time gaining eccentricity than the prograde 
orbit. Conversely, the prograde orbit spend much longer in the eccentricity decreasing phase. 

Fig. 1 illustrates the positions of the horizon, ISCO and the critical radius for prograde and retrograde orbits 
around black holes of all spins. Fig. 2 illustrates the behaviour of c for Schwarzschild orbits (a = 0) and for prograde 
and retrograde orbits around a Kerr black hole with a — .9M. The dramatic plunge in c towards negative values as 
the ISCO approaches is seen in all three cases. 



In order to calculate how the eccentricity changes as the orbit evolves, one can integrate Eq. (104) and define a new 
parameter 7, such that 



7 = ln(^) = / -dro, (111) 

ei Jr, ro 

where is the initial eccentricity at radius r^, and is the eccentricity at a smaller radius rf. Employing the numerical 
results for c, along with the analytic approximation close to the ISCO given in Eq. ( |106| ), we can numerically integrate 
this equation to derive 7. Fig. 3 shows 7 for the three cases of Fig. 2, illustrating how the eccentricity changes as 
the orbit inspirals. One can see that, in the case of a black hole with spin parameter a = 0.9M, a retrograde orbit 
{q = —0.9) will have an order-of-magnitude greater eccentricity (relative to the eccentricity it had at ro = lOOM) 
when it reaches the critical radius (the turning point on the curve) then will a prograde orbit when it reaches its 
critical radius, much further in. The amount the eccentricity increases by after the critical radius is passed depends 
crucially on the details of the physical size and mass of the orbiting particle, which it is beyond the scope of this 
paper to analyse. For a test particle with vanishing fi/M the eccentricity increases arbitrarily, but in a physical case 
this process will be cut off by the onset of dynamical instability at some point. 

X. CONCLUSIONS 

The results of this paper broadly confirm the experience of the non-rotating case, in that radiation reaction tends to 
reduce orbital eccentricity until near the the ISCO, when the onset of dynamical instability is prefigured by a period 
of decircularization of the inspiralling orbit. It seems reasonable to suppose that this effect is induced by alterations 
in the shape of the radial potential R as the ISCO approaches, since at the ISCO, the minimum which defines the 
particle's circular orbit disappears. Beyond this point the particle can only plunge towards the central body and is 
not longer in a dynamically stable orbit. One can imagine that as this point approaches, the potential well in which 
the orbit sits becomes shallower and broader (as it turns into a saddle point), so that the orbital eccentricity increases 
despite the circularizing force which drives the orbit towards the potential minimum. The tendancy of prograde orbits 
around rapidly rotating black holes to begin increasing in eccentricity only very shortly before the plunge into the 
black hole (at risco) suggests that massive bodies in such orbits will have smaller eccentricities at the end of their 
inspiral than with bodies in retrograde orbits, or the non-rotating case. In the case of prograde orbits around an 
extreme Kerr black hole, the fact that c is positive arbitrarily close to r = M, suggests that the critical radius has 
descended in the "throat" of the black hole along with the ISCO, a region where the Boyer-Lindquist co-ordinates 
become degenerate. Since our notion of circularity is so dependant on this co-ordinate system, it is unclear whether 
we can attach any meaning to the critical radius for nearly circular orbits in this extreme context. Nevertheless, from 
a practical point of view this critical radius continues to be distinguishable, in terms of the B-L radius, from the 
radius of the ISCO as we approach arbitrarily close to the case of extremal rotation, albeit that it approaches the 
latter more and more closely as the rotation increases (for prograde orbits). It is worth empasizing that our definition 
of the eccentricity, although closely tied to a particular co-ordinate system, is nevertheless an important observable 



c 



element of the gravitational wave signal emitted by the system, as seen above in Eq. (102). 

Another effect of the back reaction force on the orbit is one which tends to alter the inclination angle, which measures 
the maximum departure of the orbit from the equatorial plane. Ryan |^ has shown that nearly equatorial prograde 
orbits tend to increase their inclination angle under radiation reaction, thus moving away from being equatorial, 
although the effect is not very pronounced. Retrograde orbits, on the other hand, tend to decrease their inclination 
angle (since the spin-orbit interaction is attractive for retrograde orbits). Therefore, by the late stages of inspiral, one 
might not expect prograde orbits to have remained very close to the equatorial plane. This illustrates the need for a 
more general calculation of orbital evolution in the Kerr geometry, which deals with the issue of the Carter constant. 
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APPENDIX 

The potential functions F[r) and U{r) of the Sasaki-Nakamura equation ( |6l| ) are given in this appendix. 

n-) = -^^ (112) 



where 



and 



where 



»7 = Co + ci/r + Crjr^ + ca/r^ + C4/r'' (113) 



Co = -VlwM + A(A + 2) - \2ai^{auj - m) (114) 

ci = 8ia[3auj — X{auj — m)] (115) 

C2 = -24iaM{aLU-m) + Ua^ll - 2{auj - mf] (116) 

C3 = 24m^(at^ - to) - 24Ma2 (117) 

C4 = 12a4. (118) 

U{r)^-^^ + G' + 4^-FG (119) 

^ ' (r2 + a2)2 r2 + a2 ^ ' 



G — Ti 7^ h -—^ TTTT^ (120) 

+ a {r^ + a ) 

U, = V+ ^[(2a + |:),, - ^(a + ^)] (121) 



_,_A + 3,x^ + A+— (122) 



KB 6A 

A r 
2A 

(3 ^ 2A{-iK + r - M ). (123) 
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Tcrit/M 


-0.9 


9.64 


-0.5 


8.37 


n n 


D.Oo 


0.5 


4.70 


0.7 


3.76 


0.9 


2.56 


0.95 


2.03 


.99 


1.47 


1.0 


1.0 



TABLE I. The position of the 
critical radius, rcrit in units of M, 
for different black hole spins a. The 
parameter q = a/M is defined here 
to be negative for retrograde orbits 
and positive for prograde orbits. 




FIG. 1. Graphs showing the positions of the horizon (r+), innermost stable circular orbit (nsco) and critical radius (fcrit) 
in terms of the mean orbital radius ro for all black hole spins (o < M). Positive a indicates a prograde orbit, and negative a 
a retrograde orbit. The figures corresponding to the squares on the critical radius curve (derived numerically) are given in the 
accompanying table. 
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FIG. 2. Curves showing the evolution of the parameter c, defined in Eq. ( |l04| ) , as the mean orbital radius ro decreases, for 
three different types of orbit. For a black hole with spin a — 0.9M, both the prograde {q = a/M = 0.9) and retrograde orbits 
{q = —0.9) are shown. Also shown is the case of a Schwarzschild black hole {q = 0). In each case c begins to fall quickly 
towards zero as the innermost stable cirular orbit approaches. 
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log To/M 

FIG. 3. Curves showing the change of the orbital eccentricity as the radius ro decreases, in Schwarzschild (g = 0), and for 
prograde (g = 0.9) and retrograde {q = —0.9) orbits around a Kerr black hole with a = 0.9M. In this graph, the parameter 7 
is the natural log of the ratio of the current eccentricity (at ro) to the eccentricty the orbit had at ro = lOOM. We can see here 
clearly that at a certain point (equivalent to the critical radius illustrated in Fig. 1), the eccentricity begins to increase. For an 
arbitrarily small mass ratio fi/M it will increase indefinitely before reaching the ISCO, but in any practical case this process 
will be cut off before too long by the onset of the dynamical instability which causes the orbiting body to plunge inwards 
toward the central black hole. 
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